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ABSTRACT 


A  program  designated  FLASH  has  been  developed  to  provide  a  means 
to  study  the  propagation  of  monochromatic  radiation  from  a  plane  parallel 
source  through  a  spherical  shell  atmosphere. 

A  simple  illustration  of  the  backward  Monte  Carlo  method  utilized 
in  the  development  of  the  FLASH  program  is  discussed  in  order  tp  show 
the  advantage  ,of  the  method. 

A  brief  description  of  the  methods  employed  in  the  FLASH  program  is 
given  to  illustrate  the  application  of  the  backward  Monte  Carlo  treatment 
of  light  propagation  through  a  spherical  shell  atmosphere. 

Several  comparisons  of  FLASH  generated  data  with  data  from  other 
calculations!  methods  are  shown. 
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I.  INTRODUCTIOH 


During  the  contract  period,  1  May  1967  to  31  January  1970,  work  was 
performed  on  five  major  work  areas.  These  areas  include:  (1)  the  develop¬ 
ment  of  the  LITE-II  Monte  Carlo  procedure  for  study  of  sunlight,  moonlight, 
and  night  light  transport  in  a  plane  parallel  atmosphere  and  the  application 
of  the  LITE-II  Monte  Carlo  procedure  to  calculations  for  two  model  atmos¬ 
pheres;  (2)  the  development  of  the  LITE-I  Monte  Carlo  procedure  for  study 
of  atmospheric  transport  of  light  emitted  fay  point  monodirectional  sources 
and  the  application  of  the  LITE-I  Monte  Carlo  procedure  to  the  analysis  of 
measured  searchlignt  scattering  data;  (3)  the  application  of  the  LITE-I 
Monte  Carlo  procedure  to  the  study  of  laser  radiation  backscattering  from 
the  ground  surface  and  the  scattering  of  laser  radiation  by  the  atmosphere; 

(4)  an  analysis  of  the  computer  calculations  performed  under  Item  1  to  give 
contrast  transmission  data  for  the  two  model  atmospheres  considered;  and 

(5)  the  development  of  the  FLASH  Monte  Carlo  procedure  to  treat  sunlight 
transport  in  a  spherical  geometry  atmosphere  and  the  application  of  the  FLASH 
procedure  to  determine  skylight  distributions  for  twilight  scattering  problems. 

The  work  performed  under  Items  1  through  4  above  have  been  previously 
reported  in  References  1  through  7.  This  report  describes  the  work  per¬ 
formed  for  Item  5. 

Most  previous  calculational  methods  developed  to  analyze  the  propagation 
of  sunlight  within  the  earth’s  atmosphere  have  utilized  plane  parallel  atmos¬ 
pheric  geometries,  Ref.  8-12.  The  plane  parallel  geometry,  however,  does 
not  adequately  represent  the  earth's  atmosphere  for  twilight  conditions. 

In  an  effort  to  provide  a  calculational  tool  to  compute  the  radiation 
intensity  at  twilight,  a  Monte  Carlo  program  has  been  developed  to  treat 
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the  propagation  of  visible  and  infrared  radiation  through  a  spherical  shell 
atmosphere.  The  program,  designated  as  FLASH,  utilizes  the  backwa:  u  Monte 
Carlo  approach  to  compute  the  scattered  light  intensity  from  sunlight  as 
a  function  of  view  angle  for  receivers  located  within  or  above  the  earth's 
atmosphere.  The  Monte  Carlo  method  was  selected  because  of  the  need  to 
evaluate  the  contrxbut ‘ons  from  multiple  scattering  and  because  of  the 
versatility  offered  in  its  application  to  complex  geometries.  The  FLASH 
program  also  takes  into  account  the  bending  of  the  sun  rays  (photon  paths) 
by  the  change  in  the  atmospheric  index  of  refraction  with  altitude. 
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lx.  MONTE  CARLO  METHOD 


In  the  application  of  th~  Monte  Carlo  method  to  oLta-^  the  solution 
to  an  integral  such  as 


f(x)  dx. 


the  common  practice  is  to  express  f(x)  as  the  product  of  two  functions, 
h(x) ,  a  density  function,  and  g(x) ,  an  estimating  function.  To  evaluate 
the  integral 

b 

f  h(x)  g (x)  dx, 


the  density  function  is  normalized  by  multiplying  by  a  constanct,  c, 
where  c  is  given  by 


J  h< 


c  I  h(x)  dx.  *  1 


A  random  value  of  x  is  sampled  from  the  normalized  density  ^unction  by 
evaluating  the  equation 


RN  = 


ch(x)  dx, 


where  a£  x^<b  and  RN  is  a  random  number  Jiawn  from  a  unirorm  distribution 
between  0  and  1.  An  estimate  of  tht  integral 


f(x)  dx 
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is  obtained  by  averaging  the  estimator  function 

EST  -  g(X;L)  /  c 

for  the  sample  values  x^.  An  average  of  g(Xj,)/c  for  a  large  number,  N, 
of  random  values  of  as  given  by 

i  y  *<*,>  /  = 

to 

is  an  estimate  of  tne  mean  value  of  the  integral. 

An  arbitrary  division  of  f(x)  into  density  and  estimating  functions 
may  result  in  an  estimating  function  that  produces  extremely  high  values 
of  g(x)/c  for  some  of  the  values  of  x  sampled  from  the  density  function. 
In  selecting  the  estimating  function,  the  objective  should  be  to  obtain 
a  function  that  remains  within  reasonable  limits  for  all  of  the  sampled 
values.  The  portion  of  tne  integrand  selected  as  the  density  function 
should  be  a  function  from  which  it  is  easy  to  obtain  sample  values.  In 
some  cases,  a  better  separation  of  the  integrand  into  density  and  esti¬ 
mating  functions  can  be  achieved  if  the  integral  expression  can  be  re¬ 
written  with  a  change  of  variables. 

2.1  Backward  Monte  Carlo 

When  using  the  Backward  Monte  Carlo  approach  tc  calculate  the  scat¬ 
tered  light  intensity  at  a  point  receiver  due  to  a  plane  parallel  source 
incident  upon  a  spherical  shell  atmosphere,  a  change  in  variables  in  the 
integral  form  of  the  transport  equation  yields  an  integrand  which  is  ex¬ 
pressed  in  terms  of  integration  variables  defined  about  the  receiver 
position  rather  than  about  the  source  position. 


It  will  be  shown  that  the  integral  expressed  in  terms  of  variables 
defined  with  respect  to  the  receiver  position  allows  the  selection  of  an 
estimating  function  that  has  less  variation  with  random  values  sampled  from 
the  density  function  than  that  which  would  be  obtained  when  the  variables 
of  integration  are  expressed  in  terms  of  the  source  position. 

A  physical  Interpretation  of  the  Monte  Carlo  solution  of  the  integral 
expressed  in  terms  of  variables  defined  about  the  receiver  is  that  photon 
histories  are  started  at  the  receiver  and  traced  backward  to  the  source 
position. 

The  backward  approach  is  particularly  advantageous  for  those  problems 
in  which  the  scattered  intensity  at  the  receiver  is  desired  for  only  a 
portion  of  the  total  solid  angle.  The  backward  approach,  in  this  case, 
will  allow  all  samples  to  be  taken  within  the  solid  angle  of  interest.  For 
a  plane  source,  the  backward  approach  provides  an  optimum  sampling  of  the 
areas  on  the  source  plane  that  have  the  greatest  possibility  of  contributing 
to  the  scattered  intensity  at  the  receiver  position. 

Care  must  be  exercised  in  the  application  of  backward  Monte  Carlo 
(starting  particle  histories  at  the  receiver  and  estimating  the  direct 
intensity  at  the  first  collision)  to  Insure  that  the  integral  equation 
solved  is  equivalent  to  the  actual  physical  situation.  The  approach  taken 
was  to  write  the  Integral  equation  which  represents  a  direct  interpretation 
of  the  single  scattered  intensity  and  then  to  change  the  variables  of 
integration  so  that  sampling  of  the  particles  random  path  may  be  started 
at  the  receiver  position.  If  the  limits  of  integration  on  the  forward 
form  of  the  scattering  integral  are  transposed  in  accordance  with  the 
change  of  variables  for  the  backward  calculation,  the  resulting  equation 
should  not  only  yield  the  same  results  as  the  original  equation,  but  it 
should  lend  itself  to  a  more  economical  solution. 


'  -  ■***“  '  .  -.'J* 
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2.2  Illustration  of  Backward  Monte  Carlo 

For  illustration,  consider  an  infinite  plane  parallel  source  embedded 
within  an  infinite  media  which  scatters  radiation  isotropically.  Suppose 
that  one  wishes  to  know  the  average  single  scattered  intensity  at  some 
distance  Zq  from  the  source  plane,  Figure  1.  The  expression  for  the  single 
scattered  intensity  may  be  written  as  follows, 


SI 


2ir  oo  oo 

J 


-E  Z 

S(p,4>)  ITe  T  (1/ 4u)— 


"V 


pdpd(jidZ 


(2.1) 


where  S(p,<f>)  is  the  source  strength  at  the  position  p,(j>  on  the  source  plans, 

£t  is  the  extinction  coefficient  of  the  media,  Z  is  the  distance  from  the 
source  plane  to  the  collision  point,  and  r  is  the  distance  from  the  collision 
to  the  receiver  position.  Equation  2.1  could  be  solved  with  the  Monte  Carlo 
method  or  some  other  method  of  numerical  integration,  but  some  decision  must 
be  made  as  to  what  the  upper  limit  of  the  integral  over  the  source  radius 
will  be  and  what  delta  p  values  will  yield  the  most  efficient  calculation. 

This  is  especially  true  of  a  Monte  Carlo  solution  of  Equation  2.1,  because 
if  one  should  select  the  position  of  the  source  photon  from  a  uniform  distri¬ 
bution  over  some  finite  portion  of  the  plane,  he  could  easily  be  concen¬ 
trating  his  sampling  at  large  radial  distances  that  produce  no  significant 
single  scattered  intensity  at  the  receiver  position. 

A  better  method  for  solving  Equation  2.1  would  be  to  change  the  variables 
of  integration  from  p  and  Z  to  0  and  r,  where  the  relationship  between  the 
variables  are  as  follows: 

p  *  r  sin© 

Z  *  Z  -  r  cos0  . 
o 
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r- 

f 


€• 
i  '■ 


f, 


i* 


t* 

i- 


i. 


t" 


i 


i 
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Since  p  and  Z  are  both  functions  of  r  and  0,  we  may  transform  Equation 
2.1  into  an  equation  expressed  in  terms  of  r  and  0  through  the  use  of  the 
following  relation: 


JJf(p,z)dpdz  -jj1 


f (p(r,0) ,Z(r,0)) 


3(r,0) 


drd0 


wnere 


3(P,Z) 


is  the  absolute  value  of  the  Jacobian. 


3(r,0) 

Solving  for  the  partial  derivative  in  the  Jacobian  we  obtain: 


ifL 

3r 


sin0 


••|q~  -1  r  cos0 


and 


3Z 

3r 

3Z 

30 


-COS0 


r  sin0 


Evaluating  the  Jacobian  yields 

Je_  _le_ 


3r 

3Z 


30 

3Z 


3r  30 

Rewriting  Equation  2.1  by  substituting  r  and  0  for  p  and  Z  gives: 


SI  = 


IT 


o 


S(p,<f>)*  e 


-E  (Z  -r  cos0)  -E  r 

(l/4ir)e  sin0d0drdij) 

(2.2) 


Equation  2.2  may  be  solved  using  the  Monte  Carlo  method  by  sampling  a 
value  of  0  from  sinod0/2,  a  value  of  6  from  dd/2v  and  a  value  of  r  from 


*'  ifA 
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E^,e  dr  and  evaluating  the  remainder  of  the  expression  under  the  integral 


S(p(r ,  6)  ,<J>)  e 


-£t(Z  -r  cos0) 
1  o 


for  the  values  sampled  for  r,0,  and  <(>. 

Note  that  the  integration  of  0  in  Equation  2,2  is  over  a  finite  range 
and  since  the  integral  of 

"V 

ZTe  dr 

has  a  value  of  one  when  integrated  from  0  to  infinity,  values  of  r  approach¬ 
ing  infinity  may  be  selected  and  thus  values  of  p  -  rsin0  are  also  unlimited. 

Note  also  that  the  density  function 


and  the  estimating  function 


-E_(Z  -rcos0) 
I  o 


both  peak  for  small  values  of  r  indicating  the  sampling  is  concentrated 
near  the  area  where  the  estimating  function  has  its  highest  value.  A 
situation  of  this  nature  produces  less  statistical  fluctuation  than  would 
be  expected  for  a  Monte  Carlo  solution  of  Equation  2.1. 
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III.  CALCULATIONS  METHODS  FOR  SCATTERED  LIGHT 


The  backward  Monte  Carlo  method  discussed  in  Section  II  has  been 
utilized  in  the  development  of  the  FLAL  s  program  to  provide  a  means  to 
study  the  transport  of  visible  and  infrared  radiation  within  a  spherical 
shell  atmosphere.  A  pure  mathematician  would  probably  be  satisfied 
with  a  presentation  of  the  methods  utilized  in  FLASH  as  a  method  of 
numerically  integrating  the  integral  form  of  the  transport  equation.  He 
would  not  necessarily  require  a  physical  intrepretation  of  the  method. 
Others,  however,  conceive  of  the  Monte  Carlo  approach  as  an  experiment 
performed  on  a  digital  computer  where  the  sampling  schemes  are  related  to 
the  probability  distributions  of  the  possible  events  that  may  occur  in 
tracing  a  quantum  of  energy  as  it  scatters  through  the  atmosphere.  Apply¬ 
ing  this  concept  to  the  solution  of  the  transport  equation  expressed  in 
terms  of  variables  defined  with  respect  to  the  receiver  positions  results 
in  tracing  the  quantum  of  energy  (photon)  in  a  backward  direction  from  the 
receiver  to  the  source  and  thus  the  term  backward  Monte  Carlo.  While  the 
concept  of  the  photon  going  backward  is  physically  unrealistic,  it  does 
provide  a  means  of  interpreting  the  detailed  calculations  that  are  per¬ 
formed  by  FLASH  when  solving  the  radiation  transport  equation. 

3.1  Polarization  Parameters 

The  initial  intensity  and  polarization  characteristics  of  the 
radiation  at  the  receiver  position  are  defined  in  terms  of  the  polariza¬ 
tion  parameters  Ij_,  It| ,  U  and  V  which  are  related  to  the  Stokes  parameter 
I,  Q,  U,  and  V  by  the  relations 

1  "  h  +  Iii 

Q  -  ij.  “  ln 
u  -  u 

V  -  V 
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These  parameters  are  defined  with  reference  to  a  plane  containing  the 
radial  through  the  receiver  position  and  the  direction  of  propagation  of 
the  light  incident  at  the  receiver.  The  sunlight  incident  upon  the  earth's 
atmosphere  is  considered  to  be  unpolarized,  therefore,  the  polarization 
parameters  are  initialized  to 

l±  -  0.5W 

III  *  0.5W 

U  »  0.0 


where  W  is  a  weight  parameter  initially  set  equal  to  a  unit  intensity 
incident  from  the  direction  of  propagation  into  the  receiver  position. 
Processes  such  as  absorption  that  affect  all  polarization  parameters  in 
the  same  manner  are  included  by  an  adjustment  to  the  weight  parameter 
while  processes  such  as  Rayleigh  and  aerosol  scattering  require  individual 
adjustment  of  the  polarization  parameters.  The  product  of  the  weight 
parameter  times  the  sum  of  the  first  two  polarization  parameters  gives 
the  photon  intensity  at  any  position  in  the  atmosphere.  The  degree  of 
polarization  at  any  position  is  given  by 


(dx-i.i)2+u2+v2)  1/2 

p  - - 

dj+V 

Since  tie  t^mospheric  parameters  which  affect  the  transport  of  light 
all  vary  with  altitude,  it  wa3  necessary  that  the  geometry  of  the  FLASH 
program  be  capable  of  providing  a  means  of  defining  the  atmospheric 
conditions  as  a  function  of  altitude. 

3.2  Problem  Geometry 

The  geometry  portion  of  the  FLASH  program  provides  for  the  division 
of  a  spherical  shell  atmosphere  into  as  many  as  100  spherical  shell 
layers  or  atmospheric  layers.  The  spherical  shell  layers  are  defined 


<*•  ■*£</**  &***&»»>«.  *&*.**.. 
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by  inputing  the  altitudes  of  the  spherical  surfaces  separating  the 
layers.  The  atmosnh* tic  conditions  at  each  boundary  are  defined  by 
inputing  the  optical  distance  from  ground  level  to  the  boundary  altitude, 
the  ratio  of  the  scattering-to-pseudo  extinction  coefficient  at  the 
boundary  altitude,  the  ratio  of  the  Rayleigh-to-scattering  coefficient 
at  the  boundary  altitude,  and  the  index  of  refraction  of  the  air  just 
above  the  boundary  altitude.  The  pseudo  extinction  coefficient  is 
defined  to  be  the  sum  of  the  Rayleigh,  aerosol,  and  Ozone  coefficients. 

If  the  wavelength  of  the  light  being  considered  is  in  the  range  in  which 
water  vapor  and  absorption  is  significant,  provisions  have  been 
made  for  inputing  amounts  of  the  accumulated  water  vapor  and  accumulated 
CO2  from  ground  level  to  the  boundary  altitude.  If  the  water  vapor  and 
CO 2  concentrations  are  required,  then  the  temperature  and  pressure  for 
each  boundary  altitude  are  also  required. 

In  addition  to  the  division  of  the  atmosphere  into  atmospheric 
layers,  the  FLASH  program  provides  for  the  division  of  the  atmosphere 
into  as  many  as  10  aerosol  scattering  regions  within  each  of  which 
the  aerosol  size  distribution  remains  constant  even  though  the  number 
density  may  vary.  An  aerosol  scattering  region  is  defined  in  the  input 
by  listing  the  altitudes  of  the  upper  and  lower  spherical  surfaces 
bounding  the  region.  A  phase  matrix  number  is  also  assigned  to  each 
aerosol  scattering  region  indicating  which  of  the  normalized  aerosol 
phase  matrices  in  the  problem  input  pertains  to  the  aerosol  size 
distribution  within  that  region. 


The  source  for  the  FLASH  program  is  assumed  to  be  an  infinite  plane 

pttittllcl  SCUi.Ce  tdu^cut  tO  tuc  Upper  dtuiOophcrS  at  the  ZeiTC  pC la. IT  pGa It xCu • 

The  receiver  positions  are  all  at  the  same  altitude  but  at  different  polar 
positions  within  the  zero  azimuthal  plane. 
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3.3  Selection  of  Receiver  View  Angles 

A  set  of  solid  angles  of  view  are  defined  in  the  problem  input 
data  by  listing  the  polar  and  azimuthal  angles  that  bound  each  solid 
angle.  The  polar  angles  are  measured  from  a  radial  through  the  receiver 
position  with  the  zero  polar  angle  along  the  outward  radial,  Fig.  2. 

The  total  number  of  histories  run  with  the  FLASH  program  are 
divided  equally  among  the  polar  angle  intervals  used  to  define  the 
receiver  solid  angles  of  view.  The  number  of  histories  run  within 
a  given  polar  angle  interval  is  divided  equally  among  the  solid  angles 
defined  by  the  azimuthal  angles  bounds  listed  for  that  polar  angle 
interval. 

It  is  not  necessary  that  view  angles  be  defined  to  cover  the  4lT 
steradians  about  the  receiver.  Sampling  may  be  confined  to  any  desired 
solid  angle  or  even  to  one  or  more  discrete  directions  with  the  method 
used  to  define  the  solid  angles  of  view. 

The  FLASH  program  is  designed  to  sample  backward  photon  histories 
from  the  solid  angles  of  view  systematical] y  and  to  printout  the 
scattered  intensity  for  each  solid  angle  of  view  before  proceeding  to 
the  next. 

At  the  start  of  each  history  before  the  history  counter  NHIST  is 
incremented,  the  index  of  the  polar  angle  interval,  LA,  and  the  azimuthal 
interval,  LAZ,  within  the  polar  angle  interval  are  determined  with  the 
following  expressions: 

LA  =  1  +  NPVA*NHIST/NHMAX 
JAZ  =  NHIST  -  (LA-1) *NHMAX/NPVA 
LAZ  **  1  +  NAZ(LA)*JAZ*NPVA/NHMAX 

where  NHMAX  is  the  total  number  of  histories  to  be  run,  NPVA  is  the 
number  of  polar  angle  intervals  used  to  define  the  solid  angles  of 
view,  and  NAZ(LA)  is  the  number  of  solid  angles  of  view  in  the  LAth 
polar  angle  interval. 
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Fig,  2.  Receiver  View  Angles 


The  discrete  directions  of  the  light  entering  the  solid  angle  of 
view  is  sampled  from  a  uniform  distribution  within  the  solid  angle  of 
view.  The  'osine,  CGAM,  of  the  polar  angle  of  the  light's  direction 
entering  the  view  angles  is  given  by 

CGAM  =  COSMI (LA) -RN* (COSMI (LA) -COSMA(LA) ) 

where  BN  is  a  random  number  selected  from  s  .uniform  distribution  between 
0.0  and  1.0,  and  COSMI(LA)  and  COSMA(LA)  are  the  cosines  of  the  lower  and 
upper  bounds  of  the  polar  angles  bounding  the  LAth  polar  angle  interval. 

The  azimuthal  direction,  ALf,  of  the  Jight  entering  the  solid  angle  of 
view  is  given  by 

ALF  ■  PHIL(LA,LAZ)  +  pw*(PHIU(LA,LAZ)  -  PHIL(LA,LAZ)) 

where  RN  is  a  random  number,  and  PHIL(LA,LAZ)  and  PHIU(LA,LAZ)  are  the 
lower  and  upper  azimuthal  angle  bounds  of  the  LAZth  solid  angles  of  view 
in  the  LAth  polar  angle  interval. 

3.4  Distance  to  Collision 

Once  the  initial  direction  of  the  photon's  backward  path  is 
determined,  the  optical  distance  through  each  of  the  atmospheric 
layers  that  the  particle  traverses  along  a  path  to  the  upper  or  lower 
surface  of  the  atmosphere  is  computed.  The  path  followed  by  the  particle 
is  bent' at  each  boundary  between  the  atmospheric  layers  due  to  the  change 
in  the  index  of  refraction  from  one  layer  to  the  next.  If  the  particle's 
path  intersects  the  ground  surface  before  reaching  the  upper  boundary  of 
the  atmosphere,  the  optical  distance  between  collisions  is  chosen  from  the 
exponential  distribution 

SRHO  -  -ln(l.O-RN)  , 

where  SRHO  is  the  sampled  path  length  and  RN  is  a  random  number  chosen 
from  a  uniform  distribut  >.n  between  0.0  and  1.0.  The  random  path  length 
is  then  compared  with  th  total  optical  distance  computed  for  the  particle's 
path  from  the  receiver  to  the  point  where  the  path  intersects  the  ground 
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surface.  If  the  sampled  pathlength  Is  greater  than  the  optical  distance 
to  the  ground  intersection,  a  reflection  is  forced  to  occur  at  the  inter¬ 
section  point.  If  the  random  pathlength  is  less  than  the  optical  distance 
to  the  ground  intersection,  a  collision  is  considered  to  occur  at  the  point 
on  the  particle’s  path  where  the  optical  distance  is  equal  to  the  sample 
pathlength. 

If  the  particle’s  path  intersects  the  upper  boundary  of  the  earth's 
atmosphere  before  intersecting  the  ground  surface,  the  random  pathlength 
between  collisions  is  chosen  from  a  truncated  exponential  distribution 
thus  forcing  a  collision  before  the  particle  escapes  the  atmosphere.  The 
random  path  is  choaen  with  the  expression 

SRHO  »  -In ( 1 , 0-RN* ( 1 - 0-exp (-ROMAX) ) ) 

where  RN  is  a  random  number  chosen  from  a  uniform  distribution  between 
0,0  and  1.0  and  ROMAX  is  the  optical  distance  along  the  particle's  path 
from  the  receiver  to  the  top  of  the  atmosphere.  When  the  collision  is 
forced  to  occur  before  the  particle  escapes  the  atmosphere,  the  particle 
weight  is  adjusted  to  remove  the  bias  introduced  by  forcing  the  collision 
with  the  expression 

WATE  =  WATE (1-exp ( -ROMAX) )  , 
where  WATE  is  the  particle  weight. 

3,5  CO2  and  Water  Vapor  Absorption 

If  the  wavelength  of  light  being  considered  is  one  for  which  carbon 
dioxide  or  water  vapor  absorption  is  significant,  the  CO  2  transmissions 
are  computed  for  a  finite  wavelength  interval  by  first  computing  the 
equivalent  COj  thickness  between  collisions  for  some  standard  temperature 
Tq  and  pressure  Pq  with  the  equation 
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where  is  the  product  of  the  CO2  density  tiroes  the  distance  traversed 
in  the  ith  atmospheric  layer,  P^  is  the  pressure  in  the  ith  atmospheric 
layer,  and  is  the  temperature  in  the  ith  atmospheric  layer.  Then, 
the  CO 2  transmission  is  computed  in  FLASH  with  the  expression  (Ref.  13-14), 

X(X)  .  e-l-83(c*k(X))0-73 

The  wavelength-dependent  parameter  k(X)  is  required  as  input  to  FLASH 

and  must  be  determined  from  some  known  data  for  a  given  P  and  T  ,  i.e. 

o  o 

Ref.  15-16. 


In  the  FLASH  program,  it  is  assumed  that  the  effect  of  Rayleigh  and 
aerosol  scattering  and  ozone  absorption  for  a  small  finite  wavelength 
interval  can  be  represented  with  monochromatic  light  for  a  wavelength 
within  the  interval,  since  the  scattering  and  absorption  coefficients 
for  these  processes  vary  fairly  smoothly  with  wavelength.  However,  since 
CC>2  absorption  varies  quite  radically  with  wavelength  even  within  small 
intervals,  up  to  20  k(A)  values  may  be  input  representing  twenty  smaller 
intervals  within  some  larger  interval  and  the  average  of  the  C02  transmission 
for  the  smaller  intervals  is  computed  to  give  the  CO2  transmissions  for  the 
larger  interval. 

FLASH  uses  a  modification  of  the  expression  developed  by  Howard  et.  al. 
(Ref.  17),  who  fitted  the  Goody  model  to  their  water  vapor  absorption  data 
to  calculate  the  water  vapor  transmission. 


T(A)  ■  exp 


2 


-1.97  Wi/Wo(X) 


i  (1  +  6.57  W1/WQ(X)  ^0  1/2j1/2 


where  T(A)  is  the  transmission  at  wavelength  X,  is  the  product  of  the 
water  vapor  density  times  the  distance  traversed  in  the  ith  atmospheric 
layer,  P.^  is  the  pressure  in  the  ith  atmospheric  layer  and  WQ(X)  is  the 
vapor  length  necessary  to  reduce  the  transmission  to  50%  at  the  pressure 
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P  and  temperature  T  .  The  water  vapor  transmission  for  some  finite  wave- 
0  o 

length  interval  is  then  computed  by  averaging  the  transmission  for  as  many 
as  20  subintervals  as  ie  the  CO2  transmission.  The  photon  weight  at  each 
collision  is  multiplied  by  the  average  CO2  and  water  vapor  transmissions 
along  the  path  between  collisions  to  account  for  the  water  vapor  and  CO2 
absorption  along  that  path. 

3.6  Direction  Before  Collision 

For  each  collision  or  reflection  that  occurs  during  the  course  of 
tracing  the  life  history  of  the  particle,  the  angle  y  between  the 
particle's  direction  before  collision  and  a  radial  through  the  collision 
point  is  computed  (Fig.  3).  The  angle  a  between  the  spherical  projection 
of  the  particle's  direction  before  collision  and  the  spherical  projection 
of  the  line  joining  the  collision  and  receiver  points  is  also  calculated. 

In  addition  to  the  direction  angles  with  respect  to  the  radial  through 
the  collision  or  reflection  point,  the  displacement  from  the  receiver 
position  is  computed  in  terms  of  the  two  position  angles  6  and  <f>.  Fig.  3 
shows  the  displacement  angles  0  and  $,  where  0  is  the  angle  between  the 
rafd'-l  through  the  receiver  and  the  radial  through  the  collision  or  reflec¬ 
tion  point  and  <}>  Is  the  angle  at  the  receiver  between  the  zero  azimuth 
and  a  spherical  projection  of  a  line  joining  the  receiver  and  the  collision 
point.  The  cosine  of  the  displacement  angle,  CTHEN,  is  given  by  the 
expression 

CTHEN  -  CTHE*')DTHE-STHE*SDTHE*CALF  , 

where  CTHE  and  STHE  are  the  cosine  and  sine  of  the  displacement  angle 
at  the  previous  collision  or  reflection,  CDTHE  and  SDTHE  are  the  cosine 
and  sir.e  of  the  angle  at  the  earth's  center  subtending  the  pathlength 
between  the  current  and  previous  collisions,  and  CALF  is  the  cosine  of 
the  angle  between  the  spherical  projection  of  the  pathlength  between 
collision  and  the  spherical  projection  of  the  line  joining  the  previous 
collision  and  the  receiver  position.  The  cosine,  CALFN,  of  the  angle 


Fig.  3.  Displacement  Angles 
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between  the  spherical  projection  of  the  pathlength  and  the  spherical  pro¬ 
jection  of  the  line  joining  the  current  collision  with  the  receiver  posi¬ 
tion  is  given  by  the  expression 

CALFN  -  (CTHE-CDTHE*CTHEN) / (SDTHE*STHEN)  , 
where  the  sine  and  cosine  values  are  defined  as  above. 

3.7  Effects  of  Changes  in  Index  of  Refraction 

The  direction  that  a  photon  travels  between  collisions  is  modified 
as  it  passes  from  one  atmospheric  layer  to  the  next  by  the  change  in  the 
index  of  refraction.  If  A  is  the  angle  between  the  photon's  direction 
and  the  radial  through  the  point  of  intersection  with  an  atmospheric  layer 
boundary,  the  angle  is  modified  by  the  expression  A'  =  sin  ’'‘((n/n' )sinA) , 
where  A'  is  the  angle  the  photon  direction  makes  with  the  radial  through 
the  intersection  point  after  crossing  the  boundary  and  n  and  n'  are  the 
indices  of  refraction  for  the  original  and  new  atmospheric  layers.  In 
the  case  that  sinA  is  approximately  1.0  and  n/n'>l  so  that  sinA'  becomes 
greater  than  1.0,  sinA'  is  set  equal  to  1.0  and  the  photon  leaves  the  spher¬ 
ical  surface  in  a  tangential  plane. 

3.8  Selection  of  the  Scattering  Event 

Each  time  a  collision  or  reflection  position  is  determined,  the  FLASH 
program  locates  the  atmospheric  layer  containing  the  collision  center  and 
interpolates  betweem  the  scattering-to-pseudo  extinction  coefficient 
ratios  and  the  Rayleigh-to-scattering  coefficient  ratios  for  the  bounds 
of  the  layer  to  find  the  scattering-to-pseudo  extinction  coefficient 
ratio  and  the  Rayleigh-to-scattering  coefficient  ratio  for  the  collision 
altitude.  The  photon  weight  is  multiplied  by  the  scattering-to-pseudo 
extinction  coefficient  ratio  to  account  for  the  fact  that  no  ozone  or 
aerosol  absorption  of  the  photon  is  allowed.  The  type  of  scattering  event 
is  determined  to  be  either  a  Rayleigh  or  an  aerosol  scattering  event  by 
comparing  a  random  number  with  the  ratio  of  the  Rayleigh-to-scattering 
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coefficient.  If  the  random  number  is  less  than  the  ratio,  the  scattering 
event  is  assumed  to  be  a  Rayleigh  event,  otherwise  it  is  assumed  to  be  an 
aerosol  scattering  event. 

3.9  Selection  of  Direction,  to  Source 

Once  the  position  of  a  collision  or  a  reflection  relative  to  a 
receiver  position  has  been  determined  and  the  direction  angles  of  the 
photon  before  collision  have  been  calculated,  then  the  direction  that  the 
photon  must  take  in  order  to  intersect  the  source  plane  at  an  angle  of  90 
degrees  is  determined.  Due  to  the  bending  of  the  photon's  path  by  the 
change  in  the  atmospheric  index  of  refraction  with  altitude,  an  iterative 
process  was  developed  to  determine  the  initial  direction  to  start  the 
photon  path  so  that  it  will  intersect  the  source  plane  at  a  right  angle. 

The  iterative  process  is  described  below. 

The  position  of  the  collision  relative  to  the  zero  polar  axis  is 
computed  with  the  expression: 

cos0  =  cos0cos0p  +  sin0sin0jjCos<|>  , 

where  0  is  the  angle  between  the  radial  through  the  collision  position  and 
the  zero  polar  axis,  0  and  <}>  are  the  polar  and  azimuthal  displacement 
angles  from  the  receiver  located  at  polar  position  ©p.  If  it  were  not  for 
the  bending  of  the  photon  path,  a  light  ray  emitted  normal  tv  the  source 
plane  so  that  it  passes  through  the  collision  center  would  intersect  the 
radial  through  the  collision  center  at  the  angle  0,  Figure  4.  The  iterative 
process  used  to  determine  the  angle  of  intersection  between  the  radial  and 
the  curved  path  from  the  source  plane  to  collision  first  starts  the  photon 
path  from  the  collision  position  toward  the  source  plane  in  the  direction 
0,  Figure  4.  The  angle  of  intersection  of  the  curved  path  with  the 
source  plane  is  then  determined.  The  difference  between  the  angle  of  inter¬ 
section  and  90  degrees,  e,  is  computed  and  this  difference  is  tested  to 
determine  whether  the  angle  of  intersection  is  sufficiently  close  to  90 
degrees  to  be  considered  an  intersection  at  right  angles.  If  e  is  too 


on  of  Scattering  Angle  between  Direction  to 
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large  it  is  subtracted  from  the  angle  g  and  a  photon  path  from  the  collision 
point  is  started  in  the  direction  of  g'=g-e.  The  angle  of  intersection  of  the 
second  photon  path  and  the  source  plane  is  computed  and  tested  as  discussed 
above.  This  process  is  repeated  until  the  angle  between  the  intersection 
of  the  photon's  path  and  the  source  plane  is  within  an  acceptable  interval 
about  90  degrees. 

3.10  Polarization  Parameters  After  Scattering 

Next  the  scattering  angle  between  the  direction  of  propagation  before 
collision  and  the  direction  of  the  curved  path  normal  to  the  source  plane 
is  computed,  Figure  4.  The  cosine  of  the  scattering  angle,  CSA,  is  given 
by  the  expression 

CSA  *»  CBETP*CGAM  +  SBETP*SGAM*CPHIP  , 

where  CBETP  and  SBETP  are  the  cosine  and  sine  of  the  angle  between  the 
radial  through  the  collision  point  and  the  direction  of  the  curved  path 
normal  to  the  source  plane,  CGAM  and  SGAM  are  the  cosine  and  sine  of 
the  angle  between  the  photon's  direction  before  collision  and  a  radial 
through  the  collision  point,  and  CPHIP  is  the  cosine  of  the  angle  between 
the  spherical  projection  of  the  photon's  directions  before  and  after 
collision. 

In  order  to  determine  the  magnitude  of  the  four  polarization  para¬ 
meters  in  a  plane  containing  the  polar  axis  and  the  direction  of  the 
photon  through  the  source  plane  the  polarization  parameters  are  rotated 
into  the  plane  of  scatter  from  the  plane  containing  the  direction  before 
collision  and  the  radial  through  the  collision  point.  The  sine  and 
cosine  of  the  rotation  angle  for  the  polarization  parameters  is  computed 
with  the  equations 

SROT  =  SBETP *SPHIP/SSA 

CROT  ■  (CBETP-CSA*CGAM) / (SSA*SBETP)  , 

where  SSA  is  the  sine  of  the  scattering  angle  and  SPHIP  is  the  sine  of 
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the  angle  between  the  spherical  projections  of  the  photon's  directions 
before  and  after  collision. 

The  polarization  parameters  1^ ,1^,11  and  V  before  rotation  are 
transformed  to 

I,,'  »  I(|  *  CROT  2  +  Ij_*SR0T2  +  0.5*U*S2ROT 
Ij.'  »  I,f  *  SROT  2  +  Ij.  *CR0T2  -  0.5*U*S2ROT 
U'  =  Ij_  *  S2R0T  -  I((  *S2R0T  +  U*C2R0T 
V'  =  V 


after  rotation  through  the  angle  whose  sine  and  cosine  are  SROT  and 
CROT.  S2R0T  and  C2R0T  are  the  sine  and  cosine  of  an  angle  twice  the 
magnitude  of  the  rotation  angle. 


If  the  Rayleigh  process  is  chosen  for  the  scattering  event,  the 
polarization  parameters  after  scattering  through  the  angle  whose  cosine 
is  CSA  are 


IJr  =  (3/16ti)  *CSA2*I((' 
i;r  *  (3/1671)*^’ 

IT  =  (3/16tt)*CSA*U' 

=  (3/1677)  *CSA*V'  , 


where  I*  ,  I.'  ,  and  U'  and  V'  are  the  polarization  parameters  after  Rayleigh 
II  r  -••r  r  r 

scattering.  If  the  scattering  event  was  an  aerosol  event,  the  polarization 
parameters  after  scattering  are  defined  by 


Via  "  Ve)*V 

Xla  =  il(0)*IJ- 

u;  »  i3(e)u'-i4(0)v' 
v'a  -  i4(e)u,+i3(e)v'  , 


where  1'^,  I|a,  U'fl  and  V’&  are  the  polarization  parameters  after  aerosol 
scattering  and  i^(0),  ^(0) ,  i^ ( ©) »  and  i^C©)  are  four  elements  of 
the  aerosol  phase  matrix  evaluated  at  the  scattering  angle  8  ■  cos  1(CSA). 

3.11  Estimate  of  the  Scattered  Intensity 

An  estimate  of  the  scattered  light  intensity  at  a  receiver  is 
obtained  by  multiplying  the  photon  weight  after  the  collision  by  the  sum 
of  the  first  two  polarization  parameters  times  the  attenuation  along 
the  path  from  the  source  plane  to  the  collision. 

EST  -  MATE*  ( I,"+Ip  *exp  (-ROMAX)  , 

where  WATE  is  the  photon  weight  after  collision,  Ij J  and  Ij.'  are  the  paral¬ 
lel  and  perpendicular  components  of  the  polarization  parameters  scattered 
into  the  direction  of  the  curved  path  that  intersects  the  source  plane  at 
right  angles,  and  ROMAX  is  the  optical  distance  along  the  path  from  the 
collision  to  the  source  plane, 

3.12  Selection  of  Direction  to  Previous  Collision 

After  the  estimate  of  th  cattered  light  intensity  at  the  receiver 
from  a  given  collision  is  mas.  ,  the  FLASH  program  selects  an  angle  of 
scatter  in  order  to  trace  the  photon  through  the  next  order  of  scattering. 
The  scattering  angle  depends  upon  whether  the  scattering  event  was 
determined  to  be  a  Rayleigh  or  an  aerosol  event. 

For  either  process,  the  orientation  of  the  scattering  plane  is 
selected  from  a  uniform  distribution  between  0  and  2tt .  The  polarization 
parameters  are  rotated  into  the  scattering  plane  using  the  expressions: 

I'll  «  I,, cos 2<j>  +  Ij_  sin^tjj  +0.5Usin2d> 

l[  =  Ijj  sin^  +  cos^t})  -0.5Usin2(j) 

U'  =>  I^sin2<|)  -  Ij|  s±n24>  +Ucos2<J> 

V'  =>  V 
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where  I(j ,  I U',  and  V'  are  the  polarization  parameters  in  the  scattering 
plane,  1^  ,  1^,  U,  and  V  are  the  polarization  parameters  in  the  reference 
plane  before  rotation  and  4>  is  the  angle  between  the  reference  plane  and 
the  scattering  plane. 

If  the  RayJeigh  process  is  chosen  for  the  scattering  event,  the 
scattering  angle  is  chosen  from  the  density  function 

P(e)  *  3/8(l+cos20)sinQd6, 

and  the  Rayleigh  phase  matrix  is  applied  with  the  following  equations 

IJr  “  2cos26I/J/(Ii5  +  lV 

xir  “  211/(1./  +  IJ.) 

=  2U'cos0/(I|j  +  Ij.) 

Vj.  -  2V' cosO/d/,  +  Ij)  . 

where  V/r5  llr’  U^j  and  are  the  polarization  parameters  after  scatter¬ 
ing  through  the  angle  0  and  I,| ,  I[,  U',  and  V’  are  the  polarization 
parameters  in  the  scattering  plane  before  scattering. 

If  the  Mie  process  is  chosen  for  the  scattering  event,  the 
scattering  angle  is  chosen  from  the  density  function 

R.i  (0)  +  R  i  (0) 

-=-= - “-= -  sin0d0, 

a 

s 

xi  X f 

where  R^=  iJ+jT  >  Rjj  =  anc*  *i(®^  anc*  ^(^  are  t*ie  e-^ements  the 

aerosol  phase  matrix  relating  to  the  perpendicular  and  parallel  components 
of  the  intensity  as  computed  with  Mie  theory.  The  aerosol  phase  matrix 
is  then  applied  to  the  polarization  parameters  through  the  following 
equations  to  yield  the  value  of  the  polarization  parameters  after 
scattering: 


T//a  =  2I||i2(e)/(i1(0)+i2(6)) 

I[a  "  2I[11(0) /(i1(0)+i2 (6)) 

u;  =  2(U ,i3(0)-V'i4(0)>/(i1(0)+i2(0)) 

«  2(U'i4(0)+V,i3(0))/(i1(0)+l2(0))  . 

The  parameters  I'  ,  i!  ,  U',  and  V’  are  defined  m  the  plane  of  scatter. 

//  3  .JL3  3  3 

The  angle  <j>'  between  the  scattering  plane  and  the  new  reference  plane 
containing  the  direction  of  propagation  after  scatter  and  the  vertical 
vector  is  calculated  with  the  trigonometric  relationships  for  spherical 
triangles.  The  four  parameters  are  then  rotated  thro*'~h  this  angle  using 
the  expressions  for  rotating  the  parameters  given  previously. 

3.13  Ground  and  Cloud  Reflection 

% 

If  the  photon  is  reflected  from  the  ground  or  a  cloud  surface,  the 
FLASH  program  selects  the  direction  after  reflection  from  a  reflection 
distribution  rather  than  from  the  Rayleigh  or  aerosol  phase  matrix. 

Once  it  has  be^n  determined  that  the  photon  intersected  a  reflection 
surface,  the  indices  of  refraction  for  the  two  materials  on  either  side 
of  the  reflection  boundary  are  checked  to  determine  whether  specular  re¬ 
flection  is  possible.  If  both  indices  are  the  same,  specular  reflection  is 
not  possible;  but  if  they  are  different,  the  reflection  coefficients  for 
the  parallel  and  perpendicular  components  of  the  specular  reflectance  are 
computed  with  the  expressions: 

r(|  =  (NIcosyR+NRcosYI)/(NIoosYR  -Nrcosyx) 

rj_  =  (-Nicosyi^Nrcosyr)/C-Nicosyi+Nrcosyr)  , 

where  N_  and  N  are  the  indices  of  refraction  for  the  materials  containing 
X  R 

the  incident  and  refracted  rays  respectively  and  yT  and  y?  are  the  angle 
of  incidence  and  the  angle  of  refraction,  respectively.  The  probability 
of  specular  reflectance  is  given  by 


2 


2 
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A  random  number  is  generated  and  if  the  random  number  is  less  than  PR, 
the  reflection  is  assumed  to  be  specular.  If  the  random  number  is  greater 
than  PR,  reflection  is  assumed  to  be  diffuse. 

If  the  reflection  is  determined  to  be  specular,  the  polarization 
parameters  are  modified  by  the  expressions: 

v.  -  riiypR 

U  -  ■dIi/PE 

U'  =  r^U/PR 
V'  ■  r^r„V/PR 

If  the  reflection  event  is  determined  not  to  be  specular,  then  the 
reflection  angle  can  be  sampled  from  an  isotropic  or  cosine  distribution 
or  from  an  arbitrary  distribution  defined  by  a  set  of  tabulated  input 
data.  Non  specular  reflected  light  is  assumed  to  be  unpolarized  so  the 
polarization  parameters  for  the  reflected  photon  are  given  by 

I,|  «*  a(Ix+  I|,)/2.0 

Ij,  =  a(Ij+  I,|)/2.0 
U'  -  0.0 
V*  -  0.0  , 

where  I^j  and  Ij_  are  the  parallel  and  perpendicular  components  of  the 
intensity  associated  with  the  incident  photon,  and  a  is  the  diffuse 
albedo  for  the  reflection  surface. 
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IV.  FLASH  CALCULATIONS 

In  order  to  validate  the  backward  Monte  Carlo  approach  used  in  the 
FLASH  program,  comparisons  of  FLASH  results  with  several  other  calculations 
have  been  made. 

4.1  Comparison  for  Rayleigh  Scattering 

The  ilrst  of  these  was  a  comparison  of  the  transmitted  intensities 
through  a  unit  optical  thickness  pure  Rayleigh  atmosphere  with  the  results 
reported  by  Coulson  et..  al.  in  Ref.  9.  For  the  FLASH  calculations 
receivers  were  positioned  at  ground  level  and  at  polar  positions  such  that 
the  uncollided  light  from  the  plane  parallel  source  would  intersect  the 
radials  through  the  receivers  at  zenith  angles  whose  cosines  were  0.92, 

0.8,  0.6,  0.4,  0.2,  and  0.1.  The  FLASH  calculated  average  scattered  intensity 
multiplied  by  2ir  is  compared  in  Table  I  with  2tt  times  the  average  scattered 
intensity  obtained  by  integrating  Coulson  et.  al.'s  transmitted  intensity 
data  over  receiver  solid  angle  for  the  same  set  of  incident  zenith  angles. 

As  one  would  expect,  the  data  for  the  two  atmospheres  are  in  reasonably 
good  agreement  until  the  incident  angle  becomes  large  and  then  the  data 
for  the  spherical  shell  geometry  tends  to  over-predict  that  for  the  plane 
geometry.  The  polar  angular  distributions  of  the  scattered  intensities  for 
the  0.0  and  0.8  ground  albedos  are  compared  for  the  incident  angles  of 

cos  ■'"0.92  and  cos  ^0.4  in  Figures  5  and  6.  Although  che  statistical  fluc¬ 
tuations  of  the  Monte  Carlo  data  are  largely  due  to  the  small  sample  size  of 
1120  histories,  the  Monte  Carlo  data  does  follow  the  same  general  trend  as 
does  Coulson' s  data. 
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Table  I 


COMPARISON  OF  SCATTERED  LIGHT  INTENSITY  TRANSMITTED  THROUGH  SPHERICAL 
AND  PLANE  RAYLEIGH  ATMOSPHERES  OF  A  UNIT  OPTICAL  THICKNESS 


Cosine  2tt  *  AVERAGE  SCATTERED  INTENSITY 

Incident  Albedo 

Angle  Spherical  Plane 


0.92 

0.80 

0.60 

0.40 

0.20 


0.0 

0.606 

0.588 

0.8 

1.411 

1.373 

0.0 

0.660 

0.641 

0.8 

1.424 

1.385 

0.0 

0.745 

0.711 

0.8 

1.426 

1.373 

0.0 

0.798 

0.736 

0.8 

1.377 

1.291 

0.0 

0.757 

0.640 

0.8 

1.213 

1.068 

0.0 

0.665 

0.535 

0.8 

1.092 

0.897 

0.10 
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COSINE  EMERGENT  ANGLE 
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4.2  Comparison  for  Aerosol  Scattering 

The  second  test  problem  to  check  the  validity  of  the  FLASH  program 
was  for  an  atmosphere  containing  aerosols  (Ref.  5).  The  Rayleigh  attenua¬ 
tion  coefficients  for  the  hazy  atmosphere  were  assumed  to  be  those  given  by 
Elterman  for  a  clear  standard  atmosphere  (Ref.  18).  Also  the  ozone  absorp¬ 
tion  coefficients  for  0.55y  light  reported  in  Ref.  5  for  the  clear  standard 
atmosphere  were  used  to  define  the  ozone  absorption  as  a  function  of 
altitude  in  the  atmosphere  described  below. 

The  aerosol  attenuation  coefficient  profile  used  for  the  hazy  atmos¬ 
phere  is  shown  in  Figure  7  for  the  0.55y  wavelength  light.  The  atmosphere 
was  defined  with  a  discontinuity  at  3  km  altitude  in  the  aerosol  attenua¬ 
tion  coefficient  profile.  Below  3  km  altitude,  the  aerosol  attenuation 
coefficient  profile  was  assumed  to  vary  exponentially  with  altitude 
according  to  the  prof  ile  reported  in  Ref.  5  after  being  normalized  to 
give  a  meteorological  range  of  3  km  at  ground  level.  Above  3  km  altitude, 
the  aerosol  attenuation  coefficient  was  assumed  to  vary  with  altitude  in 
the  manner  described  by  the  average  of  105  profiles  determined  by  Elterman 
from  searchlight  measurements  taken  in  New  Mexico  during  1963  and  1964 
(Ref.  19).  The  profile  obtained  from  the  average  of  105  profiles  was 
normalized  so  that  its  magnitude  when  extrapolated  by  use  of  the  profile 
in  Ref.  5  to  ground  level  would  give  a  ground  level  meteorological 
range  of  25  km.  The  profiles  for  0-55y  wavelength  light  above  and  below 
3  km  altitude  were  normalized  to  ground  level  meteorological  ranges  with 
use  of  the  expression 

Oa  =  (3.9/V)-ok 

where  o  is  the  aerosol  attenuation  coefficient,  V  is  the  meteorological 
a 

range  in  kilometers,  and  0^  is  the  Rayleigh  attenuation  coefficient  at 
ground  level  for  0.55y  wavelength  light. 

The  aerosol  size  distribution  below  the  haze  layer  was  taken  to  be 
2  5 

proportional  to  r  '  while  the  distribution  above  the  haze  layer  was 


35 


-3  5 

taken  to  be  proportional  to  r  '  where  r  is  the  radius  of  the  aerosol 

particle.  Aerosol  phase  matrices  were  computed  for  the  aerosol  size 

-25  -35 

distributions  proportional  to  r  *  and  r  '  by  integrating  the  phase 
matrices  calculated  for  spherical  particles  with  a  real  index  of  re¬ 
fraction  of  1.33  over  the  size  range  from  0.03  to  lOy. 

A  Lambert  type  ground  surface  was  assumed  and  scattered  intensities 
were  calculated  for  ground  albedo  values  of  0.0  and  0.8. 

The  scattered  light  intensities  computed  for  the  10  kilometer  altitude 
in  the  spherical  atmosphere  discussed  above  are  compared  in  Table  II  with 
the  scattered  intensities  computed  for  the  10  kilometer  altitude  in  a 
similar  planetary  atmosphere  with  the  LITE-IV  Monte  Carlo  program.  The 
agreement  between  the  two  sets  of  data  is  within  20  percent  for  nearly 
all  incident  angles  except  for  the  light  incident  at  87.5  degrees.  This 
is  thought  to  be  a  reasonable  comparison  considering  the  fact  that  the 
data  shown  for  the  FLASH  program  were  obtained  with  a  sample  size  of 
only  960  histories. 

4.3  Intensity  Emerging  from  Earth's  Atmosphere 

As  a  further  check  out  of  the  calculational  methods  used  in  FLASH, 
two  problems  were  run  to  calculate  the  scattered  light  intensity  at  a 
satellite  orbiting  at  300  nautical  miles  above  the  earth.  The  atmosphere, 
which  was  modeled  by  a  clear  maritime  atmosphere,  has  a  ground  level 
visibility  of  25  kilometers.  The  variation  of  the  aerosol  number  density 
with  altitude  was  taken  to  be  that  given  by  the  Elterman  1968  model.  The 
aerosol  size  distribution  was  assumed  to  be  Deirmendjian's  Haze  M  model. 
Calculations  were  made  for  satellite  positions  defined  in  terms  of  the 
earth  angle  which  is  the  angle  between  the  radials  through  the  satellite 
position  and  the  sun's  position.  The  intensities  for  0.37  micron  and 
0.78  micron  wavelength  light  when  the  ground  albedo  was  taken  to  be  0.0 
and  0.8  are  compared  with  the  results  from  the  RRA-89  program  (Ref.  21) 
in  Figures  8  and  9. 
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Table  II 


COMPARISONS  OF  SCATTERED  LIGHT  AT  10  KILOMETER  ALTITUDE  IN 

ATMOSPHERE  FROM  PLANE  AND  SPHERICAL  GEOMETRY  CALCULATIONS 

A  HAZY 

Incident 

Angle 

Albedo 

2n  *  AVERAGE  SCATTERED  INTENSITY 

(Degrees) 

Spherical 

Plane 

0 

0.0 

0.411 

0.345 

0.8 

1.498 

1.533 

15 

0.0 

0.419 

0.350 

0.8 

1.459 

1.575 

30 

0.0 

0.504 

0.432 

0.8 

1.500 

1.606 

45 

0.0 

0.806 

0.605 

0.8 

1.769 

1.684 

60 

0.0 

1.009 

1.008 

0.8 

1.848 

1.922 

75 

0.0 

1.848 

1.808 

0.8 

2.333 

2.443 

87.5 

0.0 

2.698 

1.900 

0.8 


2.990 


2.184 
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The  RRA-89  procedure  was  originally  developed  for  the  purpose  of 
providing  a  method  of  computing  the  monochromatic  light  intensities  above 
model  atmospheres  that  result  from  atmospheric  scattering  and  diffuse 
reflection  by  the  ground  surface.  The  calculational  method  used  in  RRA-89 
is  based  on  the  use  of  atmospheric  reflection  data  generated  by  use  of  the 
LITE-II  procedure  for  plane  parallel  atmospheres  in  a  machine  procedure 
that  determines  the  polar  and  azimuthal  angle  distributions  of  the  light 
intensities  at  a  satellite  located  above  the  atmosphere  for  arbitrary 
satellite  position  and  direction  to  the  sun.  The  comparisons  shown  in 
Figures  8  and  9  indicate  that  results  obtained  from  the  FLASH  procedure 
are  in  good  agreement  with  data  from  RRA-89  calculations. 

4.4  Twilight  Studies 

The  FLASH  program  was  utilized  in  a  study  of  twilight  radiation  for 
wavelengths  of  0.5  and  0.7  microns  for  three  atmospheric  models.  The 
three  model  atmospheres  all  contain  the  Rayleigh  and  ozone  profiles  given 
by  Elterman  in  Reference  18.  The  aerosol  profiles  for  the  model  atmos¬ 
pheres  are  shown  in  Fig.  10  for  0,55  micron  wavelength  light.  Profile  A 
is  the  Elterman.  1964  aerosol  profile  for  a  clear  atmosphere.  Profile  B 
is  a  modification  of  the  Elterman  1968  profile  for  normal  volcanic  dust, 
and  profile  C  is  a  strong  enhancement  of  the  volcanic  dust  in  the  10  to 
20  kilometers  altitude  range  of  the  Elterman  1968  profile.  The  three  pro¬ 
files  were  normalized  to  the  ground  level  aerosol  coefficient  for  0.5p  and 
0,7).;  given  by  Elteri.-.in  in  Ref.  18.  The  scattered  light  intensity  was 
calculated  for  collimated  receivers  at  ten  meters  altitude  and  pointing  20 
degrees  elevation  in  the  plane  of  the  sun.  The  receivers  were  placed  at 
sun  depression  angles  of  0,  2,  4,  5,  6,  and  7  degrees.  Figures  11  and  12 
show  the  single  and  multiple  scattered  intensity  as  a  function  of  sun 
depression  angle  for  the  0,7  micron  wavelength  light  in  Atmospheres  A  and  C. 
The  difference  between  the  single  and  multiple  scattering  is  seen  to  oe 
greater  for  the  clear  atmosphere,  Fig.  11,  than  in  the  atmosphere  with  the 
high  concentrations  of  volcanic  dust.  Fig.  12, 


It  is  seen  in  Figs.  11  and  12  that  the  single  scattered  intensities 
at  the  receivers  for  Atmospheres  A  and  C  do  not  differ  significantly  for 
a  given  sun  depression  angle.  The  importance  of  increased  aerosol  scat¬ 
tering  that  takes  place  in  the  dust  layer  between  5  and  25  km  altitude  in 
Atmosphere  C  is  off  set  by  the  increased  attenuation  on  the  photon's 
path  through  the  atmosphere..  To  better  see  what  is  happening  to  the  light 
as  it  undergoes  single  collision  events  in  the  atmosphere,  the  FLASH 
procedure  was  modified  so  as  to  print  out  the  intensity  scattered  per  km 
of  altitude  as  a  function  of  altitude.  Fig. 13  shows  the  results  obtained 
for  0.7y  wavelength  light  undergoing  single  scattering  in  Atmospheres  A 
and  C.  It  is  seen  that  the  lower  attenuation  ir>  Atmosphere  A  results  in 
a  larger  single  scattered  intensity  over  that  c!:  Served  for  Atmosphere  C 
in  the  first  20  km  of  altitude.  The  importance  of  the  dust  layer  between 
15  and  25  km  altitude  in  Atmosphere  C  is  noted  by  the  large  single  scattered 
intensity  coming  from  those  altitudes.  Above  30  km  altitude  the  increased 
attentuation  of  the  single  scattered  radiation  is  noted  for  Atmosphere  C 
over  that  computed  for  Atmosphere  A  as  the  single  scattered  light  passes 
through  lower  altitudes  on  its  way  to  the  receiver.  Thus,  it  is  observed 
that  although  increased  aerosol  single  scattering  takes  place  within  the 
dust  layer  centered  around  20  km,  altitude  in  Atmosphere  C,  the  increased 
attenuation  resulting  from  that  dust  layer  significantly  reduces  the 
magnitude  of  the  single  scattered  intensity  reaching  the  receiver. 

An  analysis  of  the  FLASH  multiple  scattering  calculations  indicated 
that  most  of  the  multiple  scattered  intensity  at  the  receiver  position 
comes  from  photons  that  undergo  multiple  scattering  in  the  first  10  km. 
altitude  in  both  Atmospheres  A  and  C.  The  fact  that  the  number  of  photons 
undergoing  single  scattering  events  in  the  first  10  km  altitude  is  greater 
for  Atmosphere  A  than  for  Atmosphere  C  (Fig.  13)  explains  why  multiple 
scattering  is  more  important  for  Atmosphere  A  than  it  is  for  Atmosphere  C. 

Additional  studies  on  the  effect  of  dust  layers  and  multiple  scattering 
on  the  light  observed  at  a  receiver  as  a  function  of  sun  depression  angle, 
receiver  view  angle,  and  wavelength  are  being  continued  under  contract 
F19628-70-C-0156 . 
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Fig.  11 .  Scattered  Light  Intensity  at  20  Degrees  Elevation  versus  Sun  Depression  Angle  for  Aerosol 
Profile  A:  0.7ft  Wavelength  Light 
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